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^ ■ Abstract 

A convenient tool to obtain numerical methods specially tuned on oscillating 
functions is exponential fitting. Such methods are needed in various branches of 
natural sciences, particularly in physics, since a lot of physical phenomena exhibit 
a pronounced oscillatory behavior. Many exponentially-fitted (EF) symmetric 
multistep methods for y" = /(x, y) are already developed. To have an idea of the 
QQ , accuracy we examine their phase properties. The remarkably simple expression of 

the phase-lag error obtained in Theorem [2] allows to draw quantitative conclusions 
_ on the merits of each EF version. 
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1 Introduction 



Significant efforts were undertaken over the years to promote linear multistep methods 

J J 

y^'+j ^^^H n = 0, 1, . . . , (1.1) 

j=0 j=Q 

where yn+j is an approximation to y{xn+j) and fn+j = f{xn+j,yn+j), as highly com- 
petitive solvers for the special second-order initial value problem 

y" = f{x,y), y{xo) = yo, y'{xo) = y'o- (1.2) 

The method (11.11) is characterized by the polynomials p and a where 

J J 

p(C) = E"^<'' ^(C) = E^^<'' 

j=0 j=0 

We associate with method (11.11) the following linear functional 

C[h,Si]z{x) = z{x + jh) - Ypjz"{x + j h), (1.3) 

j=0 j=0 

where a is the vector of the coefficients a = [ao, . . . , aj, Po, . . . , Its algebraic order is 
defined to be p and its error constant to be Cp+2 if, for an adequately smooth arbitrary 
test function z{x), 

L[h,ai]z{x) = Cp+2hP+^z^P+^\x) + 0{hP+^). (1.4) 

In particular, 

j=0 j=0 ^' 3=0 3=0 

The principal local truncation error (plte) is the leading term of (11.41) . i.e., 

plte = Cp+2hP+^z^P+^\x). 

Throughout, we shall assume that the method satisfies the following hypotheses: 

J 

I. aj = l, \ao\ + m^O, El^J-I^O- 

3=0 
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II. p and a have no common factors. 

III. p(l) = p'(l) = and p"(l) = 2(t(1); this is necessary and sufficient for the 
method to be consistent, that is, to have order at least one. 

IV. The method is zero-stable; that is, all the roots of p lie in or on the unit circle, 
those on the unit circle having multiplicity not greater than one. 

The method is then convergent, see Henrici [1], and the polynomial p has root of 
multiplicity precisely two at +1. 

Algorithm (II. ip is said to be symmetric when 

aj = aj^j, Pj = Pj-j for all j. 

The algebraic order p and the step number J are then even numbers [2]. Symmet- 
ric multistep methods are able to preserve the amplitude of the harmonic oscillator 
y" = —uj'^y, see Lambert and Watson [2]. For the Schrodinger equation, symmetric 
two- and four-step methods received particular attention in this context [3]-[9]. For 
computations of planetary orbits with symmetric multistep methods an excellent long- 
time behavior is reported in the literature (TUl [H]. A complete explanation of the 
behavior of classical symmetric multistep methods is given in [12]. 

Usually, the coefficients of a pth-order linear multistep method are found from 
the requirement that it integrates exactly powers up to degree p+1, or equivalently, the 
operator fll.3p is vanishing for these power functions. For problems having oscillatory 
solutions, more efficient methods are obtained when they are exact for every linear 
combination of functions from the reference set 

{1, x, . . . , x^, exp(±/ia;), . . . , x'^ exp(±/ix)}, K + 2P=p — 1. (1.5) 

This technique is known as exponential fitting and has a long history [131 El- The 
set (11. 5p is characterized by two integer parameters, K and P. The set in which there 
is no classical component is identified by = — 1 while the set in which there is no 
exponential fitting component (the classical case) is identified by P = — 1. Parameter 
P will be called the level of tuning. One should take in mind that exponential fitting 
can be applied only when a good estimate of the dominant frequency of the solution 
is known in advance. The coefficients of exponentially- fitted (in short: EF) methods 
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depend on the product of the frequency /i and the stepsize h. An important property of 
EF algorithms is that they tend to the classical ones when the involved frequencies tend 
to zero, a fact which allows to say that exponential fitting represents a natural extension 
of the classical polynomial fitting. Remark that hypotheses I-IV to be convergent are 
only applicable to classical multistep methods. The examination of the convergence 
of EF multistep methods is included in Lyche's theory [T^. Many EF symmetric 
multistep methods are already constructed [3]-[n], [H]. Also, exponential fitting has 
been applied many times to other standard algorithms such as Runge-Kutta methods, 
hybrid methods, quadrature, interpolation, . . . This vast material is collected by Ixaru 
and Vanden Berghe |15j . 

To have an idea of the accuracy of the method when solving oscillatory problems 
it is more appropriate to consider the phase-lag, rather than its usual plte. We mention 
the pioneering paper of Brusa and Nigro |;16j in which the phase-lag property was 
introduced. This is actually another type of a truncation error, i.e. the angle between 
the analytical solution and the numerical solution. The purpose of this Letter is to 
investigate the phase properties of EF symmetric multistep methods. It turns out that 
for equations similar to the harmonic oscillator, the most efficient EF methods are 
those with the highest tuning level. In the case of the Schrodinger equation, this result 
was already obtained for particular two- and four-step EF multistep methods based on 
an expensive error analysis, see |6l [3, |17] . 

2 Phase-lag analysis of classical symmetric multi- 
step methods 

Linear stability analysis and phase-lag analysis of numerical methods for (11.21) is based 
on the homogeneous test equation 

y" = -uo^y, (2.6) 

where is a real constant, which may be assumed non- negative for notational con- 
venience in latter inequalities. When we apply a symmetric multistep method to the 
scalar test equation (12.61) we obtain the difference equation 

J/2 

^ji.^'^) iVn+j + Vn-j) + Ao{iy^) Vn = 0, 
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where v = uh and Aj{v'^) are polynomials in v'^. For a linear algorithm is Aj{v'^) = 
aj + i^"^ Pj, j = 0, . . . , J/2. Multistep methods with more than one stage give rise to 
higher order polynomials in z/^, see for example [12] (Sec. 5.1.6). 

The characteristic equation is 

J/2 

J2Mu'){C + Cn + Ao{i^') = 0. (2.7) 
i=i 

Of particular interest for periodic motion is the situation where those roots are on the 
unit circle. A symmetric multistep method has an interval of periodicity (0, Uq) if, for 
all u G (0, z/q), the roots Q of the characteristic equation (12.71) satisfy 

Ci = exp(a(z/)), C2 = exp(-a(z/)), |0| < 1, 3 < j < J, (2.8) 

where A(i/) is a real function of u. For any method corresponding to the characteristic 
equation (12. 7p and for all z/^ G (0, Uq), the phase-lag is defined as the difference 

t = u- X{u). 

The phase-lag order is q if 

t = cz/''+i + C>(z/''+3), (2.9) 
where c is the phase-lag constant. 

3 Phase-lag analysis of exponentially-fitted symmet- 
ric multistep methods 

Next, we consider EF symmetric multistep methods. In what follows, we consider the 
methods in their trigonometric form; i.e. instead of (11.51) . the methods are exact for 
trigonometric functions 

{1, X, ... , , cos{kx), sin(A;x), . . . cos{k x) , x^ sin(fcx)}. 

This is accomplished by setting ^ = ik m the coefficients of the EF methods. The 
coefficients are further denoted as Q.j{d) and (3j{d) where 9 = kh. The corresponding 
classical method is determined by 0:^(0) and /3j(0). An application of such an EF 
method to the scalar test equation (12.61) leads to the difference equation 

J/2 

Ajii?; 6) {yn+j + yn~j) + Ao{iy^; 6) y^ = 0. 
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The characteristic equation is 

J/2 

J2 0) iC + + Aoiiy'; 6) = 0, (3.10) 
i=i 

where Aj{v'^\ 9) are polynomials in z/^ with coefficients which depend on the parameter 
9 which specifies the method of concern. For a linear algorithm is 

A,(z/2; 9) = a,{9) + z/2 p^i^e), J = 1, . . . , J/2. 

The important modification is that the periodicity interval becomes now a two-dimensional 
region. Following Definition 6 of Coleman and Ixaru , the region of stability is a re- 
gion in the v — 9 plane, throughout which the roots of the characteristic equation (13.101) 
satisfy the periodicity condition (12.81) . As mentioned in [TJ], to study the phase prop- 
erties it is more suited to replace the pair z/, 9 by the pair i/, r = 9 /v. For any method 
corresponding to the characteristic equation (13.101) and for all (z/, r v) belonging to the 
stability region, the phase-lag is defined as (12. 9p . The phase- lag order is q if 

t = c(r)z/^+^ + 0(z/''+=^). (3.11) 

The following theorem was previously found by Simos and Williams [T8] for classical 
methods. Here, it is extended to the EF case. The proof is essentially the same as that 
given in |18] . 



Theorem 1 For all [u, r v) belonging to the stability region, a symmetric EF J -method 
with characteristic equation /l!^.l(J\) has phase-lag order q if and only if 



2E-=iA(^';^^) cos(jV) + Ao( 



z/ ; r z/ 



__c(r)z/«+2^0(z/«+4). 



Our aim is to determine the expression of c(r) for EF symmetric linear multistep 
methods. 

Firstly, we reveal a relation between the plte and the phase-lag. The analytical 
solution of (12. 6p is 

y{x) = Ci exp{iu!x) + C2 exp(— icux). 
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Some algebraic manipulation gives 

J/2 

C%a\y{x) = (2 ^A,(z/2;rz/) cos{3 u) + A^{u^-r u)) y{x). (3.12) 
i=i 

Assuming that the phase-lag order is g, i.e. (13.111) . the following identity holds (see 
Simos and Williams [18|) 

cos(j v) = cos(j A(z/)) - f c(r) u"^^ + C(z/«+^). (3.13) 

The roots of the characteristic equation (I3.10p have to satisfy the periodicity condi- 
tions ( 12. Sp . The first two conditions of ( 12. Sp are equivalent to 

J/2 

2 Aj (z/2; r z/) cos(j A(z/)) + Ao{u^; ru) = 0. (3.14) 
i=i 

We denote by plte* as the plte when solving (IM!). Using (1312|) . (l37[3l) and (EUD we 
arrive to 

J/2 

p^te* = -2c(r)z/«+2 ^ j2 ^^.^^o) y(x). (3.15) 

i=i 

A convergent symmetric method has at least order two, thus ^^=1 Q^j(O) = 2 ^j=o'^i(0)- 
Note that (t(1) = X]j=o'^i(0)- -^^ hypotheses II-III we have that cr(l) 7^ 0, so 
Sj=oi^'^i(0) 7^ 0- With this in mind, it follows from (13.150 that the phase-lag or- 
der is equal to the order of the numerical solution of (12. 6p . a result which was already 
obtained by Coleman [20] for classical two-step hybrid methods. 

Secondly, the plte of the methods considered may be written as (see [T5] ) 

plte = Cp+2 h^^^ {D^ + kY^^ y{x), (3.16) 

where = /dx^ and Cp+2 is the error constant of the corresponding classical method 
(i.e. 6' = 0). More specifically, when solving (12. 6p one easily verifies that 

plte* = (-l)f/2+i y"^^ (1 - ry+^ y{x). (3.17) 

From (I3.16p - (l3.17p it is clear that the order of the numerical solution of (12. 6p is equal to 
the order the method. Altogether, we conclude that for symmetric EF linear multistep 
methods, the phase-lag order is equal to the algebraic order. Very important to notice 
is that the algebraic order of an EF multistep method and its classical companion have 
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both the same algebraic order [131 [HI [15] . Consequently, the EF procedure conserves 
the phase-lag order. We compare fl3.15p and (13.171) to get 

2 Ei=iJ^«i(0) 

Altogether, we summarize our main observations as 

Theorem 2 An EF symmetric linear multistep method for ( fi.^) and the corresponding 
classical method have both the same phase-lag order. Furthermore, we have that 

c(r) = {l-rY^^c, 

with r = Qjv and c is the phase-lag constant of the classical method. 

When an acceptable estimate of the dominant frequency is available (i.e. r ^ 1) the 
magnitude of the phase-lag is then much smaller than that of the corresponding classical 
method (i.e. r = 0). Furthermore, the more accurate the estimate of the dominant 
frequency is, the smaller the phase-lag is. For equations similar to (12. 6p . such as the 
Schrodinger equation, it turns out that the most appropriate EF methods are those 
with the highest possible value of P. Similar results are found in [SJ [TJ |T7| for some 
two- and four-step EF methods for the Schrodinger equation via a different approach. 

As an example, we select a four-step method of Simos [9J (p. 351, Case II) 
which is determined by the parameters K = —1 and P = 3. Its order is six. The error 
constant Cg of a classical symmetric four-step method reads 

Cs = (256 + ai(0) - 3584/?2(0) - 56(3,(0) 

For this method is 

37 29 17 

ao(0) = 0, ai(0) = -l, /3o(0) = -, ^^^^^ = ^' ^^^^^ = ^' 

which results, with the help of (13.181) . in 



c(r) = (1 — r^) 

^ ^ 120960 ^ ^ 



2\4 



This can be easily checked with Theorem [H 
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